home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / hyperg_U.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  44.2 KB  |  1,407 lines
  1. /* specfunc/hyperg_U.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28. #include <gsl/gsl_sf_laguerre.h>
  29. #include <gsl/gsl_sf_pow_int.h>
  30. #include <gsl/gsl_sf_hyperg.h>
  31.  
  32. #include "error.h"
  33. #include "hyperg.h"
  34.  
  35. #define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON)
  36.  
  37. #define SERIES_EVAL_OK(a,b,x) ((fabs(a) < 5 && b < 5 && x < 2.0) || (fabs(a) <  10 && b < 10 && x < 1.0))
  38.  
  39. #define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x))
  40.  
  41.  
  42. /* Log[U(a,2a,x)]
  43.  * [Abramowitz+stegun, 13.6.21]
  44.  * Assumes x > 0, a > 1/2.
  45.  */
  46. static
  47. int
  48. hyperg_lnU_beq2a(const double a, const double x, gsl_sf_result * result)
  49. {
  50.   const double lx = log(x);
  51.   const double nu = a - 0.5;
  52.   const double lnpre = 0.5*(x - M_LNPI) - nu*lx;
  53.   gsl_sf_result lnK;
  54.   gsl_sf_bessel_lnKnu_e(nu, 0.5*x, &lnK);
  55.   result->val  = lnpre + lnK.val;
  56.   result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + 0.5*M_LNPI + fabs(nu*lx));
  57.   result->err += lnK.err;
  58.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  59.   return GSL_SUCCESS;
  60. }
  61.  
  62.  
  63. /* Evaluate u_{N+1}/u_N by Steed's continued fraction method.
  64.  *
  65.  * u_N := Gamma[a+N]/Gamma[a] U(a + N, b, x)
  66.  *
  67.  * u_{N+1}/u_N = (a+N) U(a+N+1,b,x)/U(a+N,b,x)
  68.  */
  69. static
  70. int
  71. hyperg_U_CF1(const double a, const double b, const int N, const double x,
  72.              double * result, int * count)
  73. {
  74.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  75.   const int maxiter = 20000;
  76.   int n = 1;
  77.   double Anm2 = 1.0;
  78.   double Bnm2 = 0.0;
  79.   double Anm1 = 0.0;
  80.   double Bnm1 = 1.0;
  81.   double a1 = -(a + N);
  82.   double b1 =  (b - 2.0*a - x - 2.0*(N+1));
  83.   double An = b1*Anm1 + a1*Anm2;
  84.   double Bn = b1*Bnm1 + a1*Bnm2;
  85.   double an, bn;
  86.   double fn = An/Bn;
  87.  
  88.   while(n < maxiter) {
  89.     double old_fn;
  90.     double del;
  91.     n++;
  92.     Anm2 = Anm1;
  93.     Bnm2 = Bnm1;
  94.     Anm1 = An;
  95.     Bnm1 = Bn;
  96.     an = -(a + N + n - b)*(a + N + n - 1.0);
  97.     bn =  (b - 2.0*a - x - 2.0*(N+n));
  98.     An = bn*Anm1 + an*Anm2;
  99.     Bn = bn*Bnm1 + an*Bnm2;
  100.     
  101.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  102.       An /= RECUR_BIG;
  103.       Bn /= RECUR_BIG;
  104.       Anm1 /= RECUR_BIG;
  105.       Bnm1 /= RECUR_BIG;
  106.       Anm2 /= RECUR_BIG;
  107.       Bnm2 /= RECUR_BIG;
  108.     }
  109.     
  110.     old_fn = fn;
  111.     fn = An/Bn;
  112.     del = old_fn/fn;
  113.     
  114.     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  115.   }
  116.   
  117.   *result = fn;
  118.   *count  = n;
  119.  
  120.   if(n == maxiter)
  121.     GSL_ERROR ("error", GSL_EMAXITER);
  122.   else
  123.     return GSL_SUCCESS;
  124. }
  125.  
  126.  
  127. /* Large x asymptotic for  x^a U(a,b,x)
  128.  * Based on SLATEC D9CHU() [W. Fullerton]
  129.  *
  130.  * Uses a rational approximation due to Luke.
  131.  * See [Luke, Algorithms for the Computation of Special Functions, p. 252]
  132.  *     [Luke, Utilitas Math. (1977)]
  133.  *
  134.  * z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
  135.  *
  136.  * This assumes that a is not a negative integer and
  137.  * that 1+a-b is not a negative integer. If one of them
  138.  * is, then the 2F0 actually terminates, the above
  139.  * relation is an equality, and the sum should be
  140.  * evaluated directly [see below].
  141.  */
  142. static
  143. int
  144. d9chu(const double a, const double b, const double x, gsl_sf_result * result)
  145. {
  146.   const double EPS   = 8.0 * GSL_DBL_EPSILON;  /* EPS = 4.0D0*D1MACH(4)   */
  147.   const int maxiter = 500;
  148.   double aa[4], bb[4];
  149.   int i;
  150.  
  151.   double bp = 1.0 + a - b;
  152.   double ab = a*bp;
  153.   double ct2 = 2.0 * (x - ab);
  154.   double sab = a + bp;
  155.   
  156.   double ct3 = sab + 1.0 + ab;
  157.   double anbn = ct3 + sab + 3.0;
  158.   double ct1 = 1.0 + 2.0*x/anbn;
  159.  
  160.   bb[0] = 1.0;
  161.   aa[0] = 1.0;
  162.  
  163.   bb[1] = 1.0 + 2.0*x/ct3;
  164.   aa[1] = 1.0 + ct2/ct3;
  165.   
  166.   bb[2] = 1.0 + 6.0*ct1*x/ct3;
  167.   aa[2] = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3;
  168.  
  169.   for(i=4; i<maxiter; i++) {
  170.     int j;
  171.     double c2;
  172.     double d1z;
  173.     double g1, g2, g3;
  174.     double x2i1 = 2*i - 3;
  175.     ct1   = x2i1/(x2i1-2.0);
  176.     anbn += x2i1 + sab;
  177.     ct2   = (x2i1 - 1.0)/anbn;
  178.     c2    = x2i1*ct2 - 1.0;
  179.     d1z   = 2.0*x2i1*x/anbn;
  180.     
  181.     ct3 = sab*ct2;
  182.     g1  = d1z + ct1*(c2+ct3);
  183.     g2  = d1z - c2;
  184.     g3  = ct1*(1.0 - ct3 - 2.0*ct2);
  185.     
  186.     bb[3] = g1*bb[2] + g2*bb[1] + g3*bb[0];
  187.     aa[3] = g1*aa[2] + g2*aa[1] + g3*aa[0];
  188.     
  189.     if(fabs(aa[3]*bb[0]-aa[0]*bb[3]) < EPS*fabs(bb[3]*bb[0])) break;
  190.     
  191.     for(j=0; j<3; j++) {
  192.       aa[j] = aa[j+1];
  193.       bb[j] = bb[j+1];
  194.     }
  195.   }
  196.   
  197.   result->val = aa[3]/bb[3];
  198.   result->err = 8.0 * GSL_DBL_EPSILON * fabs(result->val);
  199.   
  200.   if(i == maxiter) {
  201.     GSL_ERROR ("error", GSL_EMAXITER);
  202.   }
  203.   else {
  204.     return GSL_SUCCESS;
  205.   }
  206. }
  207.  
  208.  
  209. /* Evaluate asymptotic for z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
  210.  * We check for termination of the 2F0 as a special case.
  211.  * Assumes x > 0.
  212.  * Also assumes a,b are not too large compared to x.
  213.  */
  214. static
  215. int
  216. hyperg_zaU_asymp(const double a, const double b, const double x, gsl_sf_result *result)
  217. {
  218.   const double ap = a;
  219.   const double bp = 1.0 + a - b;
  220.   const double rintap = floor(ap + 0.5);
  221.   const double rintbp = floor(bp + 0.5);
  222.   const int ap_neg_int = ( ap < 0.0 && fabs(ap - rintap) < INT_THRESHOLD );
  223.   const int bp_neg_int = ( bp < 0.0 && fabs(bp - rintbp) < INT_THRESHOLD );
  224.  
  225.   if(ap_neg_int || bp_neg_int) {
  226.     /* Evaluate 2F0 polynomial.
  227.      */
  228.     double mxi = -1.0/x;
  229.     double nmax = -(int)(GSL_MIN(ap,bp) - 0.1);
  230.     double tn  = 1.0;
  231.     double sum = 1.0;
  232.     double n   = 1.0;
  233.     double sum_err = 0.0;
  234.     while(n <= nmax) {
  235.       double apn = (ap+n-1.0);
  236.       double bpn = (bp+n-1.0);
  237.       tn  *= ((apn/n)*mxi)*bpn;
  238.       sum += tn;
  239.       sum_err += 2.0 * GSL_DBL_EPSILON * fabs(tn);
  240.       n += 1.0;
  241.     }
  242.     result->val  = sum;
  243.     result->err  = sum_err;
  244.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(nmax)+1.0) * fabs(sum);
  245.     return GSL_SUCCESS;
  246.   }
  247.   else {
  248.     return d9chu(a,b,x,result);
  249.   }
  250. }
  251.  
  252.  
  253. /* Evaluate finite sum which appears below.
  254.  */
  255. static
  256. int
  257. hyperg_U_finite_sum(int N, double a, double b, double x, double xeps,
  258.                     gsl_sf_result * result)
  259. {
  260.   int i;
  261.   double sum_val;
  262.   double sum_err;
  263.  
  264.   if(N <= 0) {
  265.     double t_val = 1.0;
  266.     double t_err = 0.0;
  267.     gsl_sf_result poch;
  268.     int stat_poch;
  269.  
  270.     sum_val = 1.0;
  271.     sum_err = 0.0;
  272.     for(i=1; i<= -N; i++) {
  273.       const double xi1  = i - 1;
  274.       const double mult = (a+xi1)*x/((b+xi1)*(xi1+1.0));
  275.       t_val *= mult;
  276.       t_err += fabs(mult) * t_err + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
  277.       sum_val += t_val;
  278.       sum_err += t_err;
  279.     }
  280.  
  281.     stat_poch = gsl_sf_poch_e(1.0+a-b, -a, &poch);
  282.  
  283.     result->val  = sum_val * poch.val;
  284.     result->err  = fabs(sum_val) * poch.err + sum_err * fabs(poch.val);
  285.     result->err += fabs(poch.val) * (fabs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val);
  286.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  287.     result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
  288.     return stat_poch;
  289.   }
  290.   else {
  291.     const int M = N - 2;
  292.     if(M < 0) {
  293.       result->val = 0.0;
  294.       result->err = 0.0;
  295.       return GSL_SUCCESS;
  296.     }
  297.     else {
  298.       gsl_sf_result gbm1;
  299.       gsl_sf_result gamr;
  300.       int stat_gbm1;
  301.       int stat_gamr;
  302.       double t_val = 1.0;
  303.       double t_err = 0.0;
  304.  
  305.       sum_val = 1.0;
  306.       sum_err = 0.0;
  307.       for(i=1; i<=M; i++) {
  308.         const double mult = (a-b+i)*x/((1.0-b+i)*i);
  309.         t_val *= mult;
  310.     t_err += t_err * fabs(mult) + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
  311.         sum_val += t_val;
  312.     sum_err += t_err;
  313.       }
  314.  
  315.       stat_gbm1 = gsl_sf_gamma_e(b-1.0, &gbm1);
  316.       stat_gamr = gsl_sf_gammainv_e(a,  &gamr);
  317.  
  318.       if(stat_gbm1 == GSL_SUCCESS) {
  319.         gsl_sf_result powx1N;
  320.     int stat_p = gsl_sf_pow_int_e(x, 1-N, &powx1N);
  321.         double pe_val = powx1N.val * xeps;
  322.     double pe_err = powx1N.err * fabs(xeps) + 2.0 * GSL_DBL_EPSILON * fabs(pe_val);
  323.         double coeff_val = gbm1.val * gamr.val * pe_val;
  324.     double coeff_err = gbm1.err * fabs(gamr.val * pe_val)
  325.                          + gamr.err * fabs(gbm1.val * pe_val)
  326.                          + fabs(gbm1.val * gamr.val) * pe_err
  327.                          + 2.0 * GSL_DBL_EPSILON * fabs(coeff_val);
  328.  
  329.     result->val  = sum_val * coeff_val;
  330.     result->err  = fabs(sum_val) * coeff_err + sum_err * fabs(coeff_val);
  331.     result->err += 2.0 * GSL_DBL_EPSILON * (M+2.0) * fabs(result->val);
  332.     result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
  333.     return stat_p;
  334.       }
  335.       else {
  336.         result->val = 0.0;
  337.     result->err = 0.0;
  338.         return stat_gbm1;
  339.       }
  340.     }
  341.   }
  342. }
  343.  
  344.  
  345. /* Based on SLATEC DCHU() [W. Fullerton]
  346.  * Assumes x > 0.
  347.  * This is just a series summation method, and
  348.  * it is not good for large a.
  349.  *
  350.  * I patched up the window for 1+a-b near zero. [GJ]
  351.  */
  352. static
  353. int
  354. hyperg_U_series(const double a, const double b, const double x, gsl_sf_result * result)
  355. {
  356.   const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */
  357.   const double SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON;
  358.  
  359.   if(fabs(1.0 + a - b) < SQRT_EPS) {
  360.     /* Original Comment: ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X
  361.      */
  362.     /* We can however do the following:
  363.      * U(a,b,x) = U(a,a+1,x) when 1+a-b=0
  364.      * and U(a,a+1,x) = x^(-a).
  365.      */
  366.     double lnr = -a * log(x);
  367.     int stat_e =  gsl_sf_exp_e(lnr, result);
  368.     result->err += 2.0 * SQRT_EPS * fabs(result->val);
  369.     return stat_e;
  370.   }
  371.   else {
  372.     double aintb = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) );
  373.     double beps  = b - aintb;
  374.     int N = aintb;
  375.     
  376.     double lnx  = log(x);
  377.     double xeps = exp(-beps*lnx);
  378.  
  379.     /* Evaluate finite sum.
  380.      */
  381.     gsl_sf_result sum;
  382.     int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum);
  383.  
  384.  
  385.     /* Evaluate infinite sum. */
  386.  
  387.     int istrt = ( N < 1 ? 1-N : 0 );
  388.     double xi = istrt;
  389.  
  390.     gsl_sf_result gamr;
  391.     gsl_sf_result powx;
  392.     int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr);
  393.     int stat_powx = gsl_sf_pow_int_e(x, istrt, &powx);
  394.     double sarg   = beps*M_PI;
  395.     double sfact  = ( sarg != 0.0 ? sarg/sin(sarg) : 1.0 );
  396.     double factor_val = sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * gamr.val * powx.val;
  397.     double factor_err = fabs(gamr.val) * powx.err + fabs(powx.val) * gamr.err
  398.                       + 2.0 * GSL_DBL_EPSILON * fabs(factor_val);
  399.  
  400.     gsl_sf_result pochai;
  401.     gsl_sf_result gamri1;
  402.     gsl_sf_result gamrni;
  403.     int stat_pochai = gsl_sf_poch_e(a, xi, &pochai);
  404.     int stat_gamri1 = gsl_sf_gammainv_e(xi + 1.0, &gamri1);
  405.     int stat_gamrni = gsl_sf_gammainv_e(aintb + xi, &gamrni);
  406.     int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);
  407.     int stat_gamall = GSL_ERROR_SELECT_4(stat_sum, stat_gam123, stat_pochai, stat_powx);
  408.  
  409.     gsl_sf_result pochaxibeps;
  410.     gsl_sf_result gamrxi1beps;
  411.     int stat_pochaxibeps = gsl_sf_poch_e(a, xi-beps, &pochaxibeps);
  412.     int stat_gamrxi1beps = gsl_sf_gammainv_e(xi + 1.0 - beps, &gamrxi1beps);
  413.  
  414.     int stat_all = GSL_ERROR_SELECT_3(stat_gamall, stat_pochaxibeps, stat_gamrxi1beps);
  415.  
  416.     double b0_val = factor_val * pochaxibeps.val * gamrni.val * gamrxi1beps.val;
  417.     double b0_err =  fabs(factor_val * pochaxibeps.val * gamrni.val) * gamrxi1beps.err
  418.                    + fabs(factor_val * pochaxibeps.val * gamrxi1beps.val) * gamrni.err
  419.            + fabs(factor_val * gamrni.val * gamrxi1beps.val) * pochaxibeps.err
  420.            + fabs(pochaxibeps.val * gamrni.val * gamrxi1beps.val) * factor_err
  421.                    + 2.0 * GSL_DBL_EPSILON * fabs(b0_val);
  422.  
  423.     if(fabs(xeps-1.0) < 0.5) {
  424.       /*
  425.        C  X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE
  426.        C  CAREFUL IN EVALUATING THE DIFFERENCES.
  427.        */
  428.       int i;
  429.       gsl_sf_result pch1ai;
  430.       gsl_sf_result pch1i;
  431.       gsl_sf_result poch1bxibeps;
  432.       int stat_pch1ai = gsl_sf_pochrel_e(a + xi, -beps, &pch1ai);
  433.       int stat_pch1i  = gsl_sf_pochrel_e(xi + 1.0 - beps, beps, &pch1i);
  434.       int stat_poch1bxibeps = gsl_sf_pochrel_e(b+xi, -beps, &poch1bxibeps);
  435.       double c0_t1_val = beps*pch1ai.val*pch1i.val;
  436.       double c0_t1_err = fabs(beps) * fabs(pch1ai.val) * pch1i.err
  437.                        + fabs(beps) * fabs(pch1i.val)  * pch1ai.err
  438.                        + 2.0 * GSL_DBL_EPSILON * fabs(c0_t1_val);
  439.       double c0_t2_val = -poch1bxibeps.val + pch1ai.val - pch1i.val + c0_t1_val;
  440.       double c0_t2_err =  poch1bxibeps.err + pch1ai.err + pch1i.err + c0_t1_err
  441.                        + 2.0 * GSL_DBL_EPSILON * fabs(c0_t2_val);
  442.       double c0_val = factor_val * pochai.val * gamrni.val * gamri1.val * c0_t2_val;
  443.       double c0_err =  fabs(factor_val * pochai.val * gamrni.val * gamri1.val) * c0_t2_err
  444.                      + fabs(factor_val * pochai.val * gamrni.val * c0_t2_val) * gamri1.err
  445.              + fabs(factor_val * pochai.val * gamri1.val * c0_t2_val) * gamrni.err
  446.              + fabs(factor_val * gamrni.val * gamri1.val * c0_t2_val) * pochai.err
  447.              + fabs(pochai.val * gamrni.val * gamri1.val * c0_t2_val) * factor_err
  448.                      + 2.0 * GSL_DBL_EPSILON * fabs(c0_val);
  449.       /*
  450.        C  XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS)
  451.        */
  452.       gsl_sf_result dexprl;
  453.       int stat_dexprl = gsl_sf_exprel_e(-beps*lnx, &dexprl);
  454.       double xeps1_val = lnx * dexprl.val;
  455.       double xeps1_err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(beps*lnx)) * fabs(dexprl.val)
  456.                        + fabs(lnx) * dexprl.err
  457.                        + 2.0 * GSL_DBL_EPSILON * fabs(xeps1_val);
  458.       double dchu_val = sum.val + c0_val + xeps1_val*b0_val;
  459.       double dchu_err = sum.err + c0_err
  460.                       + fabs(xeps1_val)*b0_err + xeps1_err * fabs(b0_val)
  461.                       + fabs(b0_val*lnx)*dexprl.err
  462.                       + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(c0_val) + fabs(xeps1_val*b0_val));
  463.       double xn = N;
  464.       double t_val;
  465.       double t_err;
  466.  
  467.       stat_all = GSL_ERROR_SELECT_5(stat_all, stat_dexprl, stat_poch1bxibeps, stat_pch1i, stat_pch1ai);
  468.  
  469.       for(i=1; i<2000; i++) {
  470.         const double xi  = istrt + i;
  471.         const double xi1 = istrt + i - 1;
  472.     const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps);
  473.     const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps));
  474.     const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi);
  475.     const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps));
  476.         b0_val *= b0_multiplier;
  477.     b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
  478.         c0_val  = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val;
  479.     c0_err  =  fabs(c0_multiplier_1) * c0_err
  480.              + fabs(c0_multiplier_2) * b0_err
  481.          + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON
  482.          + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON;
  483.         t_val  = c0_val + xeps1_val*b0_val;
  484.     t_err  = c0_err + fabs(xeps1_val)*b0_err;
  485.     t_err += fabs(b0_val*lnx) * dexprl.err;
  486.     t_err += fabs(b0_val)*xeps1_err;
  487.     dchu_val += t_val;
  488.     dchu_err += t_err;
  489.         if(fabs(t_val) < EPS*fabs(dchu_val)) break;
  490.       }
  491.  
  492.       result->val  = dchu_val;
  493.       result->err  = 2.0 * dchu_err;
  494.       result->err += 2.0 * fabs(t_val);
  495.       result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
  496.       result->err *= 2.0; /* FIXME: fudge factor */
  497.  
  498.       if(i >= 2000) {
  499.         GSL_ERROR ("error", GSL_EMAXITER);
  500.       }
  501.       else {
  502.         return stat_all;
  503.       }
  504.     }
  505.     else {
  506.       /*
  507.        C  X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE
  508.        C  STRAIGHTFORWARD FORMULATION IS STABLE.
  509.        */
  510.       int i;
  511.       double dchu_val;
  512.       double dchu_err;
  513.       double t_val;
  514.       double t_err;
  515.       gsl_sf_result dgamrbxi;
  516.       int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);
  517.       double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps;
  518.       double a0_err =  fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err
  519.                      + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err
  520.              + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err
  521.              + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err
  522.                      + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);
  523.       stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi);
  524.  
  525.       b0_val = xeps * b0_val / beps;
  526.       b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);
  527.       dchu_val = sum.val + a0_val - b0_val;
  528.       dchu_err = sum.err + a0_err + b0_err
  529.         + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));
  530.  
  531.       for(i=1; i<2000; i++) {
  532.         double xi = istrt + i;
  533.         double xi1 = istrt + i - 1;
  534.     double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);
  535.     double b0_multiplier = (a+xi1-beps)*x/((aintb+xi1)*(xi-beps));
  536.         a0_val *= a0_multiplier;
  537.     a0_err += fabs(a0_multiplier) * a0_err;
  538.         b0_val *= b0_multiplier;
  539.     b0_err += fabs(b0_multiplier) * b0_err;
  540.         t_val = a0_val - b0_val;
  541.     t_err = a0_err + b0_err;
  542.         dchu_val += t_val;
  543.     dchu_err += t_err;
  544.         if(fabs(t_val) < EPS*fabs(dchu_val)) break;
  545.       }
  546.  
  547.       result->val  = dchu_val;
  548.       result->err  = 2.0 * dchu_err;
  549.       result->err += 2.0 * fabs(t_val);
  550.       result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
  551.       result->err *= 2.0; /* FIXME: fudge factor */
  552.  
  553.       if(i >= 2000) {
  554.         GSL_ERROR ("error", GSL_EMAXITER);
  555.       }
  556.       else {
  557.         return stat_all;
  558.       }
  559.     }
  560.   }
  561. }
  562.  
  563.  
  564. /* Assumes b > 0 and x > 0.
  565.  */
  566. static
  567. int
  568. hyperg_U_small_ab(const double a, const double b, const double x, gsl_sf_result * result)
  569. {
  570.   if(a == -1.0) {
  571.     /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x
  572.      */
  573.     result->val  = -b + x;
  574.     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
  575.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  576.     return GSL_SUCCESS;
  577.   }
  578.   else if(a == 0.0) {
  579.     /* U(0,c+1,x) = Laguerre[c,0,x] = 1
  580.      */
  581.     result->val = 1.0;
  582.     result->err = 0.0;
  583.     return GSL_SUCCESS;
  584.   }
  585.   else if(ASYMP_EVAL_OK(a,b,x)) {
  586.     double p = pow(x, -a);
  587.     gsl_sf_result asymp;
  588.     int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
  589.     result->val  = asymp.val * p;
  590.     result->err  = asymp.err * p;
  591.     result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p;
  592.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  593.     return stat_asymp;
  594.   }
  595.   else {
  596.     return hyperg_U_series(a, b, x, result);
  597.   }
  598. }
  599.  
  600.  
  601. /* Assumes b > 0 and x > 0.
  602.  */
  603. static
  604. int
  605. hyperg_U_small_a_bgt0(const double a, const double b, const double x,
  606.                       gsl_sf_result * result,
  607.               double * ln_multiplier
  608.               )
  609. {
  610.   if(a == 0.0) {
  611.     result->val = 1.0;
  612.     result->err = 1.0;
  613.     *ln_multiplier = 0.0;
  614.     return GSL_SUCCESS;
  615.   }
  616.   else if(   (b > 5000.0 && x < 0.90 * fabs(b))
  617.           || (b >  500.0 && x < 0.50 * fabs(b))
  618.     ) {
  619.     int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier);
  620.     if(stat == GSL_EOVRFLW)
  621.       return GSL_SUCCESS;
  622.     else
  623.       return stat;
  624.   }
  625.   else if(b > 15.0) {
  626.     /* Recurse up from b near 1.
  627.      */
  628.     double eps = b - floor(b);
  629.     double b0  = 1.0 + eps;
  630.     gsl_sf_result r_Ubm1;
  631.     gsl_sf_result r_Ub;
  632.     int stat_0 = hyperg_U_small_ab(a, b0,     x, &r_Ubm1);
  633.     int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub);
  634.     double Ubm1 = r_Ubm1.val;
  635.     double Ub   = r_Ub.val;
  636.     double Ubp1;
  637.     double bp;
  638.  
  639.     for(bp = b0+1.0; bp<b-0.1; bp += 1.0) {
  640.       Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
  641.       Ubm1 = Ub;
  642.       Ub   = Ubp1;
  643.     }
  644.     result->val  = Ub;
  645.     result->err  = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub);
  646.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub);
  647.     *ln_multiplier = 0.0;
  648.     return GSL_ERROR_SELECT_2(stat_0, stat_1);
  649.   }
  650.   else {
  651.     *ln_multiplier = 0.0;
  652.     return hyperg_U_small_ab(a, b, x, result);
  653.   }
  654. }
  655.  
  656.  
  657. /* We use this to keep track of large
  658.  * dynamic ranges in the recursions.
  659.  * This can be important because sometimes
  660.  * we want to calculate a very large and
  661.  * a very small number and the answer is
  662.  * the product, of order 1. This happens,
  663.  * for instance, when we apply a Kummer
  664.  * transform to make b positive and
  665.  * both x and b are large.
  666.  */
  667. #define RESCALE_2(u0,u1,factor,count)      \
  668. do {                                       \
  669.   double au0 = fabs(u0);                   \
  670.   if(au0 > factor) {                       \
  671.     u0 /= factor;                          \
  672.     u1 /= factor;                          \
  673.     count++;                               \
  674.   }                                        \
  675.   else if(au0 < 1.0/factor) {              \
  676.     u0 *= factor;                          \
  677.     u1 *= factor;                          \
  678.     count--;                               \
  679.   }                                        \
  680. } while (0)
  681.  
  682.  
  683. /* Specialization to b >= 1, for integer parameters.
  684.  * Assumes x > 0.
  685.  */
  686. static
  687. int
  688. hyperg_U_int_bge1(const int a, const int b, const double x,
  689.                   gsl_sf_result_e10 * result)
  690. {
  691.   if(a == 0) {
  692.     result->val = 1.0;
  693.     result->err = 0.0;
  694.     result->e10 = 0;
  695.     return GSL_SUCCESS;
  696.   }
  697.   else if(a == -1) {
  698.     result->val  = -b + x;
  699.     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
  700.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  701.     result->e10  = 0;
  702.     return GSL_SUCCESS;
  703.   }
  704.   else if(b == a + 1) {
  705.     /* U(a,a+1,x) = x^(-a)
  706.      */
  707.     return gsl_sf_exp_e10_e(-a*log(x), result);
  708.   }
  709.   else if(ASYMP_EVAL_OK(a,b,x)) {
  710.     const double ln_pre_val = -a*log(x);
  711.     const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
  712.     gsl_sf_result asymp;
  713.     int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
  714.     int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
  715.                                               asymp.val, asymp.err,
  716.                                               result);
  717.     return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
  718.   }
  719.   else if(SERIES_EVAL_OK(a,b,x)) {
  720.     gsl_sf_result ser;
  721.     const int stat_ser = hyperg_U_series(a, b, x, &ser);
  722.     result->val = ser.val;
  723.     result->err = ser.err;
  724.     result->e10 = 0;
  725.     return stat_ser;
  726.   }
  727.   else if(a < 0) {
  728.     /* Recurse backward from a = -1,0.
  729.      */
  730.     int scale_count = 0;
  731.     const double scale_factor = GSL_SQRT_DBL_MAX;
  732.     gsl_sf_result lnm;
  733.     gsl_sf_result y;
  734.     double lnscale;
  735.     double Uap1 = 1.0;     /* U(0,b,x)  */
  736.     double Ua   = -b + x;  /* U(-1,b,x) */
  737.     double Uam1;
  738.     int ap;
  739.  
  740.     for(ap=-1; ap>a; ap--) {
  741.       Uam1 = ap*(b-ap-1.0)*Uap1 + (x+2.0*ap-b)*Ua;
  742.       Uap1 = Ua;
  743.       Ua   = Uam1;
  744.       RESCALE_2(Ua,Uap1,scale_factor,scale_count);
  745.     }
  746.  
  747.     lnscale = log(scale_factor);
  748.     lnm.val = scale_count*lnscale;
  749.     lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val);
  750.     y.val = Ua;
  751.     y.err = 4.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(Ua);
  752.     return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  753.   }
  754.   else if(b >= 2.0*a + x) {
  755.     /* Recurse forward from a = 0,1.
  756.      */
  757.     int scale_count = 0;
  758.     const double scale_factor = GSL_SQRT_DBL_MAX;
  759.     gsl_sf_result r_Ua;
  760.     gsl_sf_result lnm;
  761.     gsl_sf_result y;
  762.     double lnscale;
  763.     double lm;
  764.     int stat_1 = hyperg_U_small_a_bgt0(1.0, b, x, &r_Ua, &lm);  /* U(1,b,x) */
  765.     int stat_e;
  766.     double Uam1 = 1.0;  /* U(0,b,x) */
  767.     double Ua   = r_Ua.val;
  768.     double Uap1;
  769.     int ap;
  770.  
  771.     Uam1 *= exp(-lm);
  772.  
  773.     for(ap=1; ap<a; ap++) {
  774.       Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
  775.       Uam1 = Ua;
  776.       Ua   = Uap1;
  777.       RESCALE_2(Ua,Uam1,scale_factor,scale_count);
  778.     }
  779.  
  780.     lnscale = log(scale_factor);
  781.     lnm.val = lm + scale_count * lnscale;
  782.     lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale));
  783.     y.val  = Ua;
  784.     y.err  = fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
  785.     y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ua);
  786.     stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  787.     return GSL_ERROR_SELECT_2(stat_e, stat_1);
  788.   }
  789.   else {
  790.     if(b <= x) {
  791.       /* Recurse backward either to the b=a+1 line
  792.        * or to a=0, whichever we hit.
  793.        */
  794.       const double scale_factor = GSL_SQRT_DBL_MAX;
  795.       int scale_count = 0;
  796.       int stat_CF1;
  797.       double ru;
  798.       int CF1_count;
  799.       int a_target;
  800.       double lnU_target;
  801.       double Ua;
  802.       double Uap1;
  803.       double Uam1;
  804.       int ap;
  805.  
  806.       if(b < a + 1) {
  807.         a_target = b-1;
  808.     lnU_target = -a_target*log(x);
  809.       }
  810.       else {
  811.         a_target = 0;
  812.     lnU_target = 0.0;
  813.       }
  814.  
  815.       stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
  816.  
  817.       Ua   = 1.0;
  818.       Uap1 = ru/a * Ua;
  819.       for(ap=a; ap>a_target; ap--) {
  820.         Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
  821.         Uap1 = Ua;
  822.         Ua   = Uam1;
  823.         RESCALE_2(Ua,Uap1,scale_factor,scale_count);
  824.       }
  825.  
  826.       if(Ua == 0.0) {
  827.         result->val = 0.0;
  828.         result->err = 0.0;
  829.     result->e10 = 0;
  830.         GSL_ERROR ("error", GSL_EZERODIV);
  831.       }
  832.       else {
  833.         double lnscl = -scale_count*log(scale_factor);
  834.         double lnpre_val = lnU_target + lnscl;
  835.         double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl));
  836.         double oUa_err   = 2.0 * (fabs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua);
  837.         int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err,
  838.                                   1.0/Ua, oUa_err,
  839.                                   result);
  840.         return GSL_ERROR_SELECT_2(stat_e, stat_CF1);
  841.       }
  842.     }
  843.     else {
  844.       /* Recurse backward to near the b=2a+x line, then
  845.        * determine normalization by either direct evaluation
  846.        * or by a forward recursion. The direct evaluation
  847.        * is needed when x is small (which is precisely
  848.        * when it is easy to do).
  849.        */
  850.       const double scale_factor = GSL_SQRT_DBL_MAX;
  851.       int scale_count_for = 0;
  852.       int scale_count_bck = 0;
  853.       int a0 = 1;
  854.       int a1 = a0 + ceil(0.5*(b-x) - a0);
  855.       double Ua1_bck_val;
  856.       double Ua1_bck_err;
  857.       double Ua1_for_val;
  858.       double Ua1_for_err;
  859.       int stat_for;
  860.       int stat_bck;
  861.       gsl_sf_result lm_for;
  862.  
  863.       {
  864.         /* Recurse back to determine U(a1,b), sans normalization.
  865.          */
  866.         double ru;
  867.     int CF1_count;
  868.         int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
  869.         double Ua   = 1.0;
  870.         double Uap1 = ru/a * Ua;
  871.         double Uam1;
  872.         int ap;
  873.         for(ap=a; ap>a1; ap--) {
  874.           Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
  875.           Uap1 = Ua;
  876.           Ua   = Uam1;
  877.       RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
  878.         }
  879.         Ua1_bck_val = Ua;
  880.     Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs(a1-a)+CF1_count+1.0) * fabs(Ua);
  881.         stat_bck = stat_CF1;
  882.       }
  883.  
  884.       if(b == 2*a1 && a1 > 1) {
  885.         /* This can happen when x is small, which is
  886.      * precisely when we need to be careful with
  887.      * this evaluation.
  888.      */
  889.     hyperg_lnU_beq2a((double)a1, x, &lm_for);
  890.     Ua1_for_val = 1.0;
  891.     Ua1_for_err = 0.0;
  892.         stat_for = GSL_SUCCESS;
  893.       }
  894.       else if(b == 2*a1 - 1 && a1 > 1) {
  895.         /* Similar to the above. Happens when x is small.
  896.      * Use
  897.      *   U(a,2a-1) = (x U(a,2a) - U(a-1,2(a-1))) / (2a - 2)
  898.      */
  899.     gsl_sf_result lnU00, lnU12;
  900.     gsl_sf_result U00, U12;
  901.     hyperg_lnU_beq2a(a1-1.0, x, &lnU00);
  902.     hyperg_lnU_beq2a(a1,     x, &lnU12);
  903.     if(lnU00.val > lnU12.val) {
  904.       lm_for.val = lnU00.val;
  905.       lm_for.err = lnU00.err;
  906.       U00.val = 1.0;
  907.       U00.err = 0.0;
  908.       gsl_sf_exp_err_e(lnU12.val - lm_for.val, lnU12.err + lm_for.err, &U12);
  909.     }
  910.     else {
  911.       lm_for.val = lnU12.val;
  912.       lm_for.err = lnU12.err;
  913.       U12.val = 1.0;
  914.       U12.err = 0.0;
  915.       gsl_sf_exp_err_e(lnU00.val - lm_for.val, lnU00.err + lm_for.err, &U00);
  916.     }
  917.     Ua1_for_val  = (x * U12.val - U00.val) / (2.0*a1 - 2.0);
  918.     Ua1_for_err  = (fabs(x)*U12.err + U00.err) / fabs(2.0*a1 - 2.0);
  919.     Ua1_for_err += 2.0 * GSL_DBL_EPSILON * fabs(Ua1_for_val);
  920.     stat_for = GSL_SUCCESS;
  921.       }
  922.       else {
  923.         /* Recurse forward to determine U(a1,b) with
  924.          * absolute normalization.
  925.          */
  926.     gsl_sf_result r_Ua;
  927.         double Uam1 = 1.0;  /* U(a0-1,b,x) = U(0,b,x) */
  928.         double Ua;
  929.         double Uap1;
  930.         int ap;
  931.     double lm_for_local;
  932.         stat_for = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_for_local); /* U(1,b,x) */
  933.     Ua = r_Ua.val;
  934.         Uam1 *= exp(-lm_for_local);
  935.     lm_for.val = lm_for_local;
  936.     lm_for.err = 0.0;
  937.  
  938.         for(ap=a0; ap<a1; ap++) {
  939.           Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
  940.           Uam1 = Ua;
  941.           Ua   = Uap1;
  942.       RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
  943.         }
  944.         Ua1_for_val  = Ua;
  945.     Ua1_for_err  = fabs(Ua) * fabs(r_Ua.err/r_Ua.val);
  946.     Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs(a1-a0)+1.0) * fabs(Ua1_for_val);
  947.       }
  948.  
  949.       /* Now do the matching to produce the final result.
  950.        */
  951.       if(Ua1_bck_val == 0.0) {
  952.         result->val = 0.0;
  953.     result->err = 0.0;
  954.     result->e10 = 0;
  955.         GSL_ERROR ("error", GSL_EZERODIV);
  956.       }
  957.       else if(Ua1_for_val == 0.0) {
  958.         /* Should never happen. */
  959.         UNDERFLOW_ERROR_E10(result);
  960.       }
  961.       else {
  962.         double lns = (scale_count_for - scale_count_bck)*log(scale_factor);
  963.     double ln_for_val = log(fabs(Ua1_for_val));
  964.     double ln_for_err = GSL_DBL_EPSILON + fabs(Ua1_for_err/Ua1_for_val);
  965.     double ln_bck_val = log(fabs(Ua1_bck_val));
  966.     double ln_bck_err = GSL_DBL_EPSILON + fabs(Ua1_bck_err/Ua1_bck_val);
  967.         double lnr_val = lm_for.val + ln_for_val - ln_bck_val + lns;
  968.     double lnr_err = lm_for.err + ln_for_err + ln_bck_err
  969.       + 2.0 * GSL_DBL_EPSILON * (fabs(lm_for.val) + fabs(ln_for_val) + fabs(ln_bck_val) + fabs(lns));
  970.     double sgn = GSL_SIGN(Ua1_for_val) * GSL_SIGN(Ua1_bck_val);
  971.         int stat_e = gsl_sf_exp_err_e10_e(lnr_val, lnr_err, result);
  972.     result->val *= sgn;
  973.         return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
  974.       }
  975.     }
  976.   }
  977. }
  978.  
  979.  
  980. /* Handle b >= 1 for generic a,b values.
  981.  */
  982. static
  983. int
  984. hyperg_U_bge1(const double a, const double b, const double x,
  985.               gsl_sf_result_e10 * result)
  986. {
  987.   const double rinta = floor(a+0.5);
  988.   const int a_neg_integer = (a < 0.0 && fabs(a - rinta) < INT_THRESHOLD);
  989.  
  990.   if(a == 0.0) {
  991.     result->val = 1.0;
  992.     result->err = 0.0;
  993.     result->e10 = 0;
  994.     return GSL_SUCCESS;
  995.   }
  996.   else if(a_neg_integer && fabs(rinta) < INT_MAX) {
  997.     /* U(-n,b,x) = (-1)^n n! Laguerre[n,b-1,x]
  998.      */
  999.     const int n = -(int)rinta;
  1000.     const double sgn = (GSL_IS_ODD(n) ? -1.0 : 1.0);
  1001.     gsl_sf_result lnfact;
  1002.     gsl_sf_result L;
  1003.     const int stat_L = gsl_sf_laguerre_n_e(n, b-1.0, x, &L);
  1004.     gsl_sf_lnfact_e(n, &lnfact);
  1005.     {
  1006.       const int stat_e = gsl_sf_exp_mult_err_e10_e(lnfact.val, lnfact.err,
  1007.                                                       sgn*L.val, L.err,
  1008.                                                       result);
  1009.       return GSL_ERROR_SELECT_2(stat_e, stat_L);
  1010.     }
  1011.   }
  1012.   else if(ASYMP_EVAL_OK(a,b,x)) {
  1013.     const double ln_pre_val = -a*log(x);
  1014.     const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
  1015.     gsl_sf_result asymp;
  1016.     int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
  1017.     int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
  1018.                                               asymp.val, asymp.err,
  1019.                                               result);
  1020.     return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
  1021.   }
  1022.   else if(fabs(a) <= 1.0) {
  1023.     gsl_sf_result rU;
  1024.     double ln_multiplier;
  1025.     int stat_U = hyperg_U_small_a_bgt0(a, b, x, &rU, &ln_multiplier);
  1026.     int stat_e = gsl_sf_exp_mult_err_e10_e(ln_multiplier, 2.0*GSL_DBL_EPSILON*fabs(ln_multiplier), rU.val, rU.err, result);
  1027.     return GSL_ERROR_SELECT_2(stat_U, stat_e);
  1028.   }
  1029.   else if(SERIES_EVAL_OK(a,b,x)) {
  1030.     gsl_sf_result ser;
  1031.     const int stat_ser = hyperg_U_series(a, b, x, &ser);
  1032.     result->val = ser.val;
  1033.     result->err = ser.err;
  1034.     result->e10 = 0;
  1035.     return stat_ser;
  1036.   }
  1037.   else if(a < 0.0) {
  1038.     /* Recurse backward on a and then upward on b.
  1039.      */
  1040.     const double scale_factor = GSL_SQRT_DBL_MAX;
  1041.     const double a0 = a - floor(a) - 1.0;
  1042.     const double b0 = b - floor(b) + 1.0;
  1043.     int scale_count = 0;
  1044.     double lm_0, lm_1;
  1045.     double lm_max;
  1046.     gsl_sf_result r_Uap1;
  1047.     gsl_sf_result r_Ua;
  1048.     int stat_0 = hyperg_U_small_a_bgt0(a0+1.0, b0, x, &r_Uap1, &lm_0);
  1049.     int stat_1 = hyperg_U_small_a_bgt0(a0,     b0, x, &r_Ua,   &lm_1);
  1050.     int stat_e;
  1051.     double Uap1 = r_Uap1.val;
  1052.     double Ua   = r_Ua.val;
  1053.     double Uam1;
  1054.     double ap;
  1055.     lm_max = GSL_MAX(lm_0, lm_1);
  1056.     Uap1 *= exp(lm_0-lm_max);
  1057.     Ua   *= exp(lm_1-lm_max);
  1058.  
  1059.     /* Downward recursion on a.
  1060.      */
  1061.     for(ap=a0; ap>a+0.1; ap -= 1.0) {
  1062.       Uam1 = ap*(b0-ap-1.0)*Uap1 + (x+2.0*ap-b0)*Ua;
  1063.       Uap1 = Ua;
  1064.       Ua   = Uam1;
  1065.       RESCALE_2(Ua,Uap1,scale_factor,scale_count);
  1066.     }
  1067.  
  1068.     if(b < 2.0) {
  1069.       /* b == b0, so no recursion necessary
  1070.        */
  1071.       const double lnscale = log(scale_factor);
  1072.       gsl_sf_result lnm;
  1073.       gsl_sf_result y;
  1074.       lnm.val = lm_max + scale_count * lnscale;
  1075.       lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + scale_count * fabs(lnscale));
  1076.       y.val  = Ua;
  1077.       y.err  = fabs(r_Uap1.err/r_Uap1.val) * fabs(Ua);
  1078.       y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
  1079.       y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
  1080.       y.err *= fabs(lm_0-lm_max) + 1.0;
  1081.       y.err *= fabs(lm_1-lm_max) + 1.0;
  1082.       stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  1083.     }
  1084.     else {
  1085.       /* Upward recursion on b.
  1086.        */
  1087.       const double err_mult = fabs(b-b0) + fabs(a-a0) + 1.0;
  1088.       const double lnscale = log(scale_factor);
  1089.       gsl_sf_result lnm;
  1090.       gsl_sf_result y;
  1091.  
  1092.       double Ubm1 = Ua;                 /* U(a,b0)   */
  1093.       double Ub   = (a*(b0-a-1.0)*Uap1 + (a+x)*Ua)/x;    /* U(a,b0+1) */
  1094.       double Ubp1;
  1095.       double bp;
  1096.       for(bp=b0+1.0; bp<b-0.1; bp += 1.0) {
  1097.         Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
  1098.         Ubm1 = Ub;
  1099.         Ub   = Ubp1;
  1100.         RESCALE_2(Ub,Ubm1,scale_factor,scale_count);
  1101.       }
  1102.  
  1103.       lnm.val = lm_max + scale_count * lnscale;
  1104.       lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
  1105.       y.val = Ub;
  1106.       y.err  = 2.0 * err_mult * fabs(r_Uap1.err/r_Uap1.val) * fabs(Ub);
  1107.       y.err += 2.0 * err_mult * fabs(r_Ua.err/r_Ua.val) * fabs(Ub);
  1108.       y.err += 2.0 * GSL_DBL_EPSILON * err_mult * fabs(Ub);
  1109.       y.err *= fabs(lm_0-lm_max) + 1.0;
  1110.       y.err *= fabs(lm_1-lm_max) + 1.0;
  1111.       stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  1112.     }
  1113.     return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  1114.   }
  1115.   else if(b >= 2*a + x) {
  1116.     /* Recurse forward from a near zero.
  1117.      * Note that we cannot cross the singularity at
  1118.      * the line b=a+1, because the only way we could
  1119.      * be in that little wedge is if a < 1. But we
  1120.      * have already dealt with the small a case.
  1121.      */
  1122.     int scale_count = 0;
  1123.     const double a0 = a - floor(a);
  1124.     const double scale_factor = GSL_SQRT_DBL_MAX;
  1125.     double lnscale;
  1126.     double lm_0, lm_1, lm_max;
  1127.     gsl_sf_result r_Uam1;
  1128.     gsl_sf_result r_Ua;
  1129.     int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
  1130.     int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
  1131.     int stat_e;
  1132.     gsl_sf_result lnm;
  1133.     gsl_sf_result y;
  1134.     double Uam1 = r_Uam1.val;
  1135.     double Ua   = r_Ua.val;
  1136.     double Uap1;
  1137.     double ap;
  1138.     lm_max = GSL_MAX(lm_0, lm_1);
  1139.     Uam1 *= exp(lm_0-lm_max);
  1140.     Ua   *= exp(lm_1-lm_max);
  1141.  
  1142.     for(ap=a0; ap<a-0.1; ap += 1.0) {
  1143.       Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
  1144.       Uam1 = Ua;
  1145.       Ua   = Uap1;
  1146.       RESCALE_2(Ua,Uam1,scale_factor,scale_count);
  1147.     }
  1148.  
  1149.     lnscale = log(scale_factor);
  1150.     lnm.val = lm_max + scale_count * lnscale;
  1151.     lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
  1152.     y.val  = Ua;
  1153.     y.err  = fabs(r_Uam1.err/r_Uam1.val) * fabs(Ua);
  1154.     y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
  1155.     y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
  1156.     y.err *= fabs(lm_0-lm_max) + 1.0;
  1157.     y.err *= fabs(lm_1-lm_max) + 1.0;
  1158.     stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  1159.     return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  1160.   }
  1161.   else {
  1162.     if(b <= x) {
  1163.       /* Recurse backward to a near zero.
  1164.        */
  1165.       const double a0 = a - floor(a);
  1166.       const double scale_factor = GSL_SQRT_DBL_MAX;
  1167.       int scale_count = 0;
  1168.       gsl_sf_result lnm;
  1169.       gsl_sf_result y;
  1170.       double lnscale;
  1171.       double lm_0;
  1172.       double Uap1;
  1173.       double Ua;
  1174.       double Uam1;
  1175.       gsl_sf_result U0;
  1176.       double ap;
  1177.       double ru;
  1178.       double r;
  1179.       int CF1_count;
  1180.       int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
  1181.       int stat_U0;
  1182.       int stat_e;
  1183.       r = ru/a;
  1184.       Ua   = GSL_SQRT_DBL_MIN;
  1185.       Uap1 = r * Ua;
  1186.       for(ap=a; ap>a0+0.1; ap -= 1.0) {
  1187.         Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
  1188.         Uap1 = Ua;
  1189.         Ua   = Uam1;
  1190.     RESCALE_2(Ua,Uap1,scale_factor,scale_count);
  1191.       }
  1192.  
  1193.       stat_U0 = hyperg_U_small_a_bgt0(a0, b, x, &U0, &lm_0);
  1194.  
  1195.       lnscale = log(scale_factor);
  1196.       lnm.val = lm_0 - scale_count * lnscale;
  1197.       lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_0) + fabs(scale_count * lnscale));
  1198.       y.val  = GSL_SQRT_DBL_MIN*(U0.val/Ua);
  1199.       y.err  = GSL_SQRT_DBL_MIN*(U0.err/fabs(Ua));
  1200.       y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a0-a) + CF1_count + 1.0) * fabs(y.val);
  1201.       stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  1202.       return GSL_ERROR_SELECT_3(stat_e, stat_U0, stat_CF1);
  1203.     }
  1204.     else {
  1205.       /* Recurse backward to near the b=2a+x line, then
  1206.        * forward from a near zero to get the normalization.
  1207.        */
  1208.       int scale_count_for = 0;
  1209.       int scale_count_bck = 0;
  1210.       const double scale_factor = GSL_SQRT_DBL_MAX;
  1211.       const double eps = a - floor(a);
  1212.       const double a0 = ( eps == 0.0 ? 1.0 : eps );
  1213.       const double a1 = a0 + ceil(0.5*(b-x) - a0);
  1214.       gsl_sf_result lnm;
  1215.       gsl_sf_result y;
  1216.       double lm_for;
  1217.       double lnscale;
  1218.       double Ua1_bck;
  1219.       double Ua1_for;
  1220.       int stat_for;
  1221.       int stat_bck;
  1222.       int stat_e;
  1223.       int CF1_count;
  1224.  
  1225.       {
  1226.         /* Recurse back to determine U(a1,b), sans normalization.
  1227.          */
  1228.         double Uap1;
  1229.         double Ua;
  1230.         double Uam1;
  1231.         double ap;
  1232.         double ru;
  1233.         double r;
  1234.         int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
  1235.         r = ru/a;
  1236.         Ua   = GSL_SQRT_DBL_MIN;
  1237.         Uap1 = r * Ua;
  1238.         for(ap=a; ap>a1+0.1; ap -= 1.0) {
  1239.           Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
  1240.           Uap1 = Ua;
  1241.           Ua   = Uam1;
  1242.       RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
  1243.         }
  1244.         Ua1_bck = Ua;
  1245.         stat_bck = stat_CF1;
  1246.       }
  1247.       {
  1248.         /* Recurse forward to determine U(a1,b) with
  1249.          * absolute normalization.
  1250.          */
  1251.     gsl_sf_result r_Uam1;
  1252.     gsl_sf_result r_Ua;
  1253.     double lm_0, lm_1;
  1254.         int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
  1255.         int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
  1256.         double Uam1 = r_Uam1.val;
  1257.         double Ua   = r_Ua.val;
  1258.         double Uap1;
  1259.         double ap;
  1260.  
  1261.     lm_for = GSL_MAX(lm_0, lm_1);
  1262.     Uam1 *= exp(lm_0 - lm_for);
  1263.     Ua   *= exp(lm_1 - lm_for);
  1264.  
  1265.         for(ap=a0; ap<a1-0.1; ap += 1.0) {
  1266.           Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
  1267.           Uam1 = Ua;
  1268.           Ua   = Uap1;
  1269.       RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
  1270.         }
  1271.         Ua1_for = Ua;
  1272.     stat_for = GSL_ERROR_SELECT_2(stat_0, stat_1);
  1273.       }
  1274.  
  1275.       lnscale = log(scale_factor);
  1276.       lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale;
  1277.       lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs(scale_count_for - scale_count_bck)*fabs(lnscale));
  1278.       y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck;
  1279.       y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val);
  1280.       stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  1281.       return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
  1282.     }
  1283.   }
  1284. }
  1285.  
  1286.  
  1287. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  1288.  
  1289.  
  1290. int
  1291. gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x,
  1292.                              gsl_sf_result_e10 * result)
  1293. {
  1294.   /* CHECK_POINTER(result) */
  1295.  
  1296.   if(x <= 0.0) {
  1297.     DOMAIN_ERROR_E10(result);
  1298.   }
  1299.   else {
  1300.     if(b >= 1) {
  1301.       return hyperg_U_int_bge1(a, b, x, result);
  1302.     }
  1303.     else {
  1304.       /* Use the reflection formula
  1305.        * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
  1306.        */
  1307.       gsl_sf_result_e10 U;
  1308.       double ln_x = log(x);
  1309.       int ap = 1 + a - b;
  1310.       int bp = 2 - b;
  1311.       int stat_e;
  1312.       int stat_U = hyperg_U_int_bge1(ap, bp, x, &U);
  1313.       double ln_pre_val = (1.0-b)*ln_x;
  1314.       double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs(b)+1.0) * fabs(ln_x);
  1315.       ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */
  1316.       stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
  1317.                                             U.val, U.err,
  1318.                                             result);
  1319.       return GSL_ERROR_SELECT_2(stat_e, stat_U);
  1320.     }
  1321.   }
  1322. }
  1323.  
  1324.  
  1325. int
  1326. gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x,
  1327.                          gsl_sf_result_e10 * result)
  1328. {
  1329.   const double rinta = floor(a + 0.5);
  1330.   const double rintb = floor(b + 0.5);
  1331.   const int a_integer = ( fabs(a - rinta) < INT_THRESHOLD );
  1332.   const int b_integer = ( fabs(b - rintb) < INT_THRESHOLD );
  1333.  
  1334.   /* CHECK_POINTER(result) */
  1335.  
  1336.   if(x <= 0.0) {
  1337.     DOMAIN_ERROR_E10(result);
  1338.   }
  1339.   else if(a == 0.0) {
  1340.     result->val = 1.0;
  1341.     result->err = 0.0;
  1342.     result->e10 = 0;
  1343.     return GSL_SUCCESS;
  1344.   }
  1345.   else if(a_integer && b_integer) {
  1346.     return gsl_sf_hyperg_U_int_e10_e(rinta, rintb, x, result);
  1347.   }
  1348.   else {
  1349.     if(b >= 1.0) {
  1350.       /* Use b >= 1 function.
  1351.        */
  1352.       return hyperg_U_bge1(a, b, x, result);
  1353.     }
  1354.     else {
  1355.       /* Use the reflection formula
  1356.        * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
  1357.        */
  1358.       const double lnx = log(x);
  1359.       const double ln_pre_val = (1.0-b)*lnx;
  1360.       const double ln_pre_err = fabs(lnx) * 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(b));
  1361.       const double ap = 1.0 + a - b;
  1362.       const double bp = 2.0 - b;
  1363.       gsl_sf_result_e10 U;
  1364.       int stat_U = hyperg_U_bge1(ap, bp, x, &U);
  1365.       int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
  1366.                                             U.val, U.err,
  1367.                                             result);
  1368.       return GSL_ERROR_SELECT_2(stat_e, stat_U);
  1369.     }
  1370.   }
  1371. }
  1372.  
  1373.  
  1374. int
  1375. gsl_sf_hyperg_U_int_e(const int a, const int b, const double x, gsl_sf_result * result)
  1376. {
  1377.   gsl_sf_result_e10 re;
  1378.   int stat_U = gsl_sf_hyperg_U_int_e10_e(a, b, x, &re);
  1379.   int stat_c = gsl_sf_result_smash_e(&re, result);
  1380.   return GSL_ERROR_SELECT_2(stat_c, stat_U);
  1381. }
  1382.  
  1383.  
  1384. int
  1385. gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result)
  1386. {
  1387.   gsl_sf_result_e10 re;
  1388.   int stat_U = gsl_sf_hyperg_U_e10_e(a, b, x, &re);
  1389.   int stat_c = gsl_sf_result_smash_e(&re, result);
  1390.   return GSL_ERROR_SELECT_2(stat_c, stat_U);
  1391. }
  1392.  
  1393.  
  1394. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  1395.  
  1396. #include "eval.h"
  1397.  
  1398. double gsl_sf_hyperg_U_int(const int a, const int b, const double x)
  1399. {
  1400.   EVAL_RESULT(gsl_sf_hyperg_U_int_e(a, b, x, &result));
  1401. }
  1402.  
  1403. double gsl_sf_hyperg_U(const double a, const double b, const double x)
  1404. {
  1405.   EVAL_RESULT(gsl_sf_hyperg_U_e(a, b, x, &result));
  1406. }
  1407.